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Abstract. 

Virus trafficking is fundamental for infection success and plasmid cytosolic trafficking is a key step 
of gene delivery. Based on the main physical properties of the cellular transport machinery such as 
microtubules, motor proteins, our goal here is to derive a mathematical model to study cytoplasmic 
trafficking. Because experimental results reveal that both active and passive movement are necessary 
for a virus to reach the cell nucleus, by taking into account the complex interactions of the virus 
with the microtubules, we derive here an estimate of the mean time a virus reaches the nucleus. In 
particular, we present a mathematical procedure in which the complex viral movement, oscillating 
between pure diffusion and a deterministic movement along microtubules, can be approximated by 
a steady state stochastic equation with a constant effective drift. An explicit expression for the drift 
amplitude is given as a function of the real drift, the density of microtubules and other physical 
parameters. The present approach can be used to model viral trafficking inside the cytoplasm, which 
is a fundamental step of viral infection, leading to viral replication and in some cases to cell damage. 
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1. Introduction. Because cytosolic transport has been identified as a critical 
barrier for synthetic gene delivery [1] , plasmids or viral DNAs delivery from the cell 
membrane to the nuclear pores has attracted the attention of many biologists. The 
cell cytosol contains many types of organelles, actin filaments, microtubules and many 
others, so that to reach the nucleus, a viral DNA has to travel through a crowded 
and risky environment. We are interested here in studying the efficiency of the de- 
livery process and we present a mathematical model of virus trafficking inside the 
cell cytoplasm. We model the viral movement as a Brownian motion. However, the 
density of actin filaments and microtubules, inside the cell, can hinder diffusion, as 
demonstrated experimentally [2]. In a crowded environment, we will model the virus 
as a material point. This reduction is simplistic for several reasons: actin filament 
network can trapped a diffusing object and beyond a certain size, as observed exper- 
imentally, a DNA fragment cannot find its way across the actin filaments [2] ■ Active 
directional transport along microtubules or actin filaments seems then the only way to 
deliver a plasmid to the nucleus. The active transport of the virus involves in general 
motor proteins, such as Kinesin (to travel in the direction of the cell membrane) or 
Dynein (to travel toward the nucleus). Once a virus is attached to a Dynein protein, 
its movement can be modeled as a deterministic drift toward the nucleus. 

Recently, a macroscopic modeling has been developed to describe the dynamics 
of adenovirus concentration inside the cell cytoplasm [3J. This approach offers very 
interesting results about the effect of microtubules, but neglects the complexity of the 
geometry and cannot be used to describe the movement of a single virus, which might 
be enough to cause cellular infection. Modeling a virus trafficking imposes to use a 
stochastic description. We model here the motion of a virus as that of a material point, 
so the probability of its trapping by actin filaments or microtubules is neglected. In 
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the present approximation, the viral movement has two main components: a Brownian 
one, which accounts for its free movement, and a drift directed towards the centrosome 
or MTOC (Microtubules Organization Center), an organelle located near the nucleus. 
The magnitude of the drift along microtubules depends on many parameters, such as 
the binding and unbinding rates and the velocity of the motor proteins [4]. 

In the present approach, we present a method to approximate a time dependent 
dynamics of virus trafficking by an effective stochastic equation with a radial steady 
state drift. The main difficulties we have to overcome arise from the time depen- 
dent nature of the trajectories which consists of intermittent epochs of drifts and free 
diffusion. We propose to derive an explicit expression for the steady state drift am- 
plitude. In this approximation, the effective drift will gather the mean properties of 
the cytoplasmic organization such as the density of microtubules and its off binding 
rate. 

Our method to find the effective drift can be described as follow: first, we ap- 
proximate the cell geometry as a two dimensional disk and use a pure Brownian 
description to approximate the virus diffusion step. This geometrical approximation 
is valid, for any two dimensional cell such as the in vitro flat skin fibroblast culture 
cells [3]: indeed, due to their adhesion to the substrate, the thickness of these cells 
can be neglected in first approximation. Second, when the distribution of the initial 
viral position is uniform on the cell surface, we will estimate, during the diffusing 
period, the hitting position on a microtubule. By solving a partial differential equa- 
tion, inside a sliced shape domain, delimited by two neighboring microtubules, we 
will provide an estimate of the mean time to the most likely hitting point. Finally, 
the amplitude of the radial steady state drift will be obtained by an iterative method 
which assumes that, after a virus has moved a certain distance along a microtubule, 
it is released at a point uniformly distributed on the final radial distance from the 
nucleus, ready for a new random walk. This scenario repeats until the virus reaches 
the nucleus surface. Finally, we will compute the mean time, the mean number of 
steps before a virus reaches the nucleus and the amplitude of the effective drift by 
using the following criteria: the Mean First Passage Time (MFPT) to the nucleus 
of the iterative approximation is equal to the MFPT obtained by solving directly an 
Ornstcin-Uhlenbeck stochastic equation. The explicit computation of the effective 
drift is a key result in the estimation of the probability and the mean time a single 
virus or DNA molecule takes to reach a small nuclear pore [5]. 

2. Modeling stochastic viral movement inside a biological cell. We ap- 
proximate the cell as a two dimensional geometrical domain f2, which is here a disk 
of radius R and the nucleus located inside is a concentric disk of much smaller radius 
S << R. We model the motion of an unattached DNA fragment as a material point, 
so that the probability of its trapping by actin filaments or microtubules is neglected. 
The motion of a (DNA) molecule of mass m is described by the overdamped limit of 
the Langevin equation (Smoluchowski's limit) [5] for the position X(f) of the molecule 
at time t. When the particle is not bound to a microtubule filament, its movement 
is described as pure Brownian with a diffusion constant D. When the particle hits 
a filament, it binds for a certain random time and moves along with a determinist 
drift. We only take into account the movement toward the nucleus, which is confound 
here with the MTOC (Microtubule organization center), an organelle where all mi- 
crotubules converge (see figure (I2.1[) ). For 5 < |X(t)| < R, we describe the overall 
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movement by the stochastic rule 

( V2Dw for X(t) free 
±={ (2.1) 
[ Vrn for X (t) bound 

where V is a constant velocity, w a <5-correlated standard white noise and r the 
X radial coordinate, the origin of which is the center of the cell. We assume that 
all filaments starting from the cell surface end on the nucleus surface. The binding 
time corresponds to a chemical reaction event and we assume that it is exponentially 
distributed and for simplicity we approximate it by a constant t m . 

Once a virus enters the cell membrane, its moves according to the rule (|2.1|) , until 
it hits a nuclear pore. Although nuclear pores occupy a small portion of the nuclear 
surface, we only consider the virus movement until it hits the nuclear surface D (5) . In 
this article, our goal is to replace equation (|2.1[) by a steady state stochastic equation 

X = b(X) + \/2£>w, (2.2) 

where the drift b is radially symmetric. In a first approximation, we consider a 
constant radial drift b(X) = — B-^ and compute hereafter the value of the constant 
amplitude B such that the MFPT of the process (|2.2[) and (|2.ip to the nucleus are 
equal. 




(a) (b) 



Fig. 2.1. Cell geometry, (a) Cell's microtubules network. All microtubules starting from the 
cell membrane converge to the Microtubule Organization center (MTOC), located near the nucleus, 
(b) simplified cell's microtubules network organization. The MTOC coincides with the nucleus. 

2.1. Modeling viral dynamics in the cytoplasm. Inside the cytosol, micro- 
tubules are distributed on the cell surface and converging radially to the MTOC. We 
denote by p this distribution (see figure (|2.ip ). We do not take into account in the 
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present analysis, the effect of organelle crowding due to the endoplasmic reticulum, 
the Golgi apparatus and many others. However, it is always possible to include them 
indirectly by using an apparent diffusion constant. We consider the fundamental 
domain fl defined as the two dimensional slice of angle O between two neighboring 
microtubules. We consider here that microtubules are uniformly distributed and thus 
= 2^, where N is the total number of microtubules. 

Although a virus can drift along microtubules in both directions by using dynein 
(resp. kinesin) motor proteins for the inward (resp. forward) movement, we only take 
into account the drift toward the nucleus [7]. It is still unclear what is the precise 
mechanism used by a virus to select a direction of motion. Attached to a dynein 
molecule, the virus transport consists in several steps of few nanometers: the length 
of each step depends on the load of the transported cargo and ATP-concentration 
[8]. We neglect here the complexity of this process, assuming that ATP molecules are 
abundant, uniformly distributed over the cell and is not a limiting factor. We thus 
assume the bound particle moves towards the nucleus with the mean constant velocity 
V. When the particle is released away from the microtubule, inside the domain, the 
process can start afresh and the particle diffuses freely. Because the Smoluchowski 
limit of the Langevin equation does not account for the change in velocity, we release 
the the particle at a certain distance away from the microtubule, but at a fixed 
distance from the nucleus (at an angle chosen uniformly distributed), see figure l2~2l 

Because microtubules are taken uniformly distributed, we can always release the 
virus inside the slice fi, between two neighboring microtubules. Thus the movement 
of the virus will be studied in Q: inside the cytosol, the viral movement is purely 
Brownian until it hits a microtubule which is now the lateral boundary of f2 (see 
figure (|2.2| ). We assume that once a virus hits a microtubule, with probability one, 
the dynamics switches from diffusion to a dctcrminist motion with a constant drift. A 
virus spends on a microtubule a time that we consider to be exponentially distributed, 
since this time is the sum of escape time from deep potential wells. We approximate 
the total time on a microtubule by the mean time t m . Thus a virus moves to a 
distance d m = Vt m along microtubule, which depends only on the characteristic of 
the virus-microtubule interactions. To summarize, the virus trajectory is a succession 
of diffusion steps mixed with some periods of attaching and detaching to microtubules. 
Thus scenario repeats until the virus hits the nucleus surface (Figure (|2.2[) ). 

2.2. Computing the MFPT to reach the nucleus. We define the mean 
time to infection as the MFPT a virus reaches the surface of the disk D (S) inside the 
domain f2 (see figure (|2.2|Q . 

To estimate the mean time to infection, we note that we can decompose the overall 
motion as a repeated fundamental step. This step consists of the free diffusion of the 
particle inside the domain followed by the motion along the microtubule. The total 
time of infection Ti is then the sum of times the particle spends in each step. Although 
the time on microtubule is determinist equal to t m , the diffusing time is not easy to 
compute and depend on the initial condition. Ultimately Ti depends on the number 
of times the fundamental step repeats before the particle reaches the nucleus. 

Let us now described each step: the first step starts when the virus enter the cell 
at the periphery r = R = Rq (at a random angle 9 G [0; 0]) and ends when the virus 
hits either the lateral boundary or the nucleus. We now consider the first passage 
time u (Rq) to the absorbing boundary and by r(Ra) the hitting position. To account 
for the determinist drift, we move during a deterministic time t m the virus from a 
distance d m along the microtubule. In that case, the initial random position for the 
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Fig. 2.2. Virus trafficking inside a cell, (a) Representation of the cell portion between 
two microtubules, (b) Transport along microtubules: Two fundamental steps are represented. A 
fundamental step is made of the two intermediate step which are first the diffusion inside the domain 
followed by the directed motion along the microtubule. 

next step is given by r = R\ = r(Ro) — d m and the total time in step 1 is U (Rq) + t m . 

We iterate the process as follow and consider in each step k the distance Rk = 
r(-Rfc_x) — d m from which the particle starts and the time u (Rk) +t rn it spends inside 
the step. If we denote by n s the random number of steps necessary to reach the 
nucleus r = S, the time to infection Tj is given by 



n»-l 



n = ^2 u{R k ) + n s t m + t r 



(2.3) 



fc=0 



where t r is a residual time, which is the time to reach the nucleus before a full step is 
completed. 

We are interested in the estimating the mean first passage MFPT t of Tj, given 

by 



t = E{n) = E ( ^ u ( R k) ) + <n s > t m + < t r >, 



(2.4) 



where < n s > is the mean number of steps and < t r > is the mean residual time. If 
we introduce the probability density function p m = Pr{n s = m} that the number of 
step is exactly equal to m, we can write 



t = E(n) = ^ E ( X] u{R k )\n s = m p m + < n s > t m + < t r >, (2.5) 



k=0 



To estimate the MFPT r, we shall approximate the previous sum by using the mean 
first passage time u(Rk) in each step k. To estimate u(Rk), we will solve (in the next 
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paragraph) the Dynkin's equation with the following boundary conditions: inside fi, 
the particle is reflected at the periphery r = R, absorbed at the nucleus dft a and at 
9 = and 9 = 0. We will also estimate the mean distance dk covered during step 
k. For that purpose we will estimate the mean exit position r m {Rk), conditioned on 
the initial position r = Rk- Indeed, we will thus get dk = Rk — r m (Rk) — d m . The 
estimates of the mean distances covered for each fundamental step will ultimately lead 
to an approximation of the mean number of step n =< n s >: n will be computed 
such that R n > 8 and R n +i < 8 (where R n = r m (i?„_i) — d m is defined recursively). 
Finally, we will obtain the following approximation for the infection time 

r w ^u{R k ) +nt m + < t r >, (2.6) 
k=o 

The mean residual time < t r > can be equal cither to u(R„) + at m , where < a < 1 
if the virus binds to a microtubule in the last step and travels a distance ad m on the 
microtubule, or to the MFPT to the nuclear boundary if r m (R n ) < 8. 

3. Mean First Passage Time and Exit point distribution . In first ap- 
proximation, under the assumptions of a sufficiently small radius 8 << R and an angle 
9 << 1 , for the computation of the MFPT and the distribution of exit points, we 
neglect the nuclear area. We define the full pie wedge tt R domain of angle O. Inside 
Sl R , we use the boundary conditions described above. Consequently, the MFPT to a 
microtubule u = u{r, 9) of a virus starting initially at position (r, 9) is solution of the 
Dynkin's equations [6] 

DAu (x) = -1 for x E n R (3.1) 
u(x) =0 for x e dfl R 

|^ = o for x e an? , 

on 

where dVL R = {9 = 0} U {9 = 6} and Q R = {r = R}. 

3.1. The general solution for the MFPT. In this paragraph only we reparametrize 
the domain by —0/2 < 9 < 0/2. By writing equation (|3.ip in polar coordinates and 
using the separation of variables, the general solution of equation 

2 u ldu 1 d 2 u\ , „ N _ „ , „ % 7} 



is given by [9] 



Or* ' rdr ' ? ^jM=-^M^ ^ 
u (r, 9) = Q for (r, 0) € dVt R . (3.3) 



u ( r > ^) — TFi { C ^T^ - 1 )+ E ^"^"cos (A„0) , for - 9/2 < < 9/2(3.4) 
4L> V cos (9 / ' 

where the edge boundary is here located at position (9 = ±9/2. The sum in the 
right-hand side is the general solution of the homogeneous problem Au = in £l R . 
The boundary conditions on the sides of the wedge impose that 

A„ = (2n+1)^, (3.5) 
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while the reflecting condition for r — R reads 
du 



(R,e) = for all e [-9/2,9/2]. (3.6) 
Or 



Using the uniqueness of Fourrier decomposition and the boundary condition (|3.6p , we 
obtain that 

_ (-l)" +1 8fl 2 - A - 

" DQXl {XI 4) • (3 ' 7) 

By averaging formula (|3.4|l over an initial uniform distribution, the MFPT to a one 
of the wedge is given by 



, , 1 . r 2 /ton (9) \ ^, 16i? 2 - A "r A " 



n=0 

where A„ = (2n + 1) For 9 small, equation (|3 . 8|) can be approximated by 



? 2 I r W© 



. (r) -^(^g-.V ^ggliil r . (3.9) 



Dtt 3 (tt/9) -4 



3.2. Exit points distribution. To estimate the position a virus will attach 
preferentially to the microtubule, we determine the distribution of exit points, when 
the viral particle initially started at a radial distance from the nucleus. We recall 
that the probability density function (pdf) p(r,t|r ) to find a diffusing particle in a 
volume element dr at time t inside the wedge f2, conditioned on the initial position 
r = ro is solution of the diffusion equation 

d p(r ^ |ro) = DAp (r, t\v ) for r e V R 
p(r,t|r ) = for r G dfl* 

o far r E , 

On 

where the initial condition is p (r, 0|ro) = 8 (r — ro)- The distribution of exit points 
e (y) is given by 

POO 

e(v)= / 3(v,t)dt, (3.10) 
J o 

where the flux j is defined by 

J [Vi t) = —D — 

On |r = y 

If we denote C (ro, r) = L°° p (r, f |ro) eft then C is solution of 

-DAC(r ,r)=<5(r-r ), (3.11) 
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and 

dC 

e (y) = -D— (r , y) for y e fi« (3.12) 

Consequently, to obtain the pdf of exit points e, we use the Green function in the wedge 
domain Vl R . By using a conformal transformation, we hereafter solve a simplified case 
of an open wedge (i.e. without a reflecting boundary at r = R). This computation 
could be compared with the general one that will be derived in the next section. 

To compute the exit points distribution, we consider the solution of equation 
(|3.11[) , obtained by the image method and a conformal transformation from the open 
wedge to the upper complex half-plane. The Green function, solution of equation (|3 . lip 
in the upper complex half-plane is given by 

C(z)^-^-ln Z —^, (3.13) 
2irL> z — Zq 

where Zq the complex conjugate of zq. Using the conformal transformation u) = 
f (z) = z& [10], that maps the interior of the wedge of opening angle to the upper 
half plane, the Green function in the wedge is given by 

C (z) = -L-ln ( Z " " Z °\ ) . (3.14) 

The flux to the line 9 is given by 

DdC , iB . 1 ip(re ie y.(k -k* ) 



r 89 v > 2Tir ((re* 6 )" - k ) ((re^f - k*) 
1 -2v{re ie ) u r%sin{v6 Q ) 



2irr 



' {re ie f v + rl u - 2 (re !9 )" r»cos {u9 ) ' 



where v — ^, fco = z o = { r o e * 9 °Y ■ Finally, the exit point distribution for 9 = is 
given by 



6 r 2v + rl v + 2 (rr ) v cos (v9 ) 



£ e M = , ^ I «w„„ v — ,..a T > ( 3 - 15 ) 



while for = it is given by 



A matlab check guarantees that 

{ee(r)+eo(r)}dr = l. (3.17) 



This simple computation is instructive and shall be compared to the full one given in 
section 13.31 
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3.3. Exit pdf in a Pie Wedge. To compute the exit points distribution in a 
pie wedge with a reflecting boundary at r — i?, we search for an explicit solution of 
the diffusion equation in polar coordinates inside the pie wedge. We first consider the 
general diffusion equation 

f(M^(^ + if + ^g§)(*,%> (3.18) 
P(x,0\y) =6(x-y) 

where the boundary conditions are given in (|3.1[) . We may often use the change of 
variable Vn e N* : 

nir 

The initial condition is given by 

2 

p (x, 0\y) = p (r, 9, 0\r , 6 ) = -tt—S (r - r Q ) } sin (kO) sin (k9 ) , 

for 9 < 6o (if > 6 , #o must be replaced by 9 — 9 ). To compute the solution of 
equation (|3.18|) . we consider the Laplace transform p of the probability p 

sp (r, 6, s\r , ) - ^—5 (r - r ) 2J sin (k9) sin (k9 ) = D + -— + ~^^J ( r ; > s l r o, ) • 

Using the separation of variables, we have 

p (r, 9, s\ro, 9q) — R k (r, s) sin (k9) sin (kOo) , 

k 

Using the change of variable, x (s) = rw^ and xq (s) = ro\fjj, we get for all k that 



1 / k 2 \ 2 

Rk i x ( s ) ! s )+— r^ R k (s) , a)- 1 + ——a Rk (x (s) , s) = -— — {x (s) - x (s)) . 

x{s) y x(s) J ei)xo(s) 

(3-19) 

Rk ( x (s) , s) is a superposition of modified Bessel functions of order k : I k (x (s)) and 
K k (x (s)) for x (s) x (s) : 

R k (x (s) , s) = A k I k {x (a)) + B k K k (x (a)) , 

where A k and B k are real constants. Since K k diverges as x (s) — > 0, the interior 
solution for (x (s) < xq (a)) depends only on I k . We denote by D k the exterior solution 
for(x (s) > xq (a)). We use the general notation x A y — min(x,y) and x V y = 
max (x, y), thus 

R k (x (s) , s) = A k I k (x (s) A x (s)) D k (x (s) V x (a)) . 

To determine D k = a k I k + b k K k , we use the reflecting condition at x (s) — x+ (s) = 
R\fl5 an< ^ we S e t that 

A k I k (x Q (a)) . (a k h. (x+ (a)) + b k K h (x+ (a))) = 0. 
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We choose 

fl fe = ~Kk ( x + ( s )) and b k = Ik ( x + (s)) • 

Thus 

R k (x (s) , s) = A k I k (x (s) A x («)) (l' k (x+ (»)) - < (z+ («)) 4) (a: (a) V x («)) . 

The constants ^4^ are determined by integrating equation (|3.19|) over an infinitesimal 
interval that includes r . Using the continuity of R k , we get 



(/4)ie(s)>xo(s) \x(s)=x (s) ~ (-^fc)x(s)<x (s) U(s)=xo(s) 



9£>a;o(s)' 
that is 

A, (l k (l' k (x + (s)) K k - K' k (x + (sj) l' k ) - 4 (l' k (x+ («)) iffc - < (x+ ( s )) J fc )) (x («)) = - 
after some simplifications, we get 

A k l' k (x + (s)) (l k K k - l' k K k ) (x (sj) = - eD l Q{s y 
Using the recurrent relation between modified Bessel functions (see [IT] or page 489 

4 ipo (s)) = ( 4-1 Tt4 ) (xo ( s )) and K k (x (s)) = ( -iffc-i TT-^fe ) (xo (s)) , 

we get 

(x+ («)) fj fc f-^ fc _x ^-r^fe) - flfc-i ^4^ #fc) (aso (*)) = -; 



ODx (s) 



x (s) / V x (s) J J ' QDx (s) ' 

that is 

2 



A k I k {x+ (a)) (/feiffe-i + Jfc_iJffc) (10 («)) 



eDx (s) ' 

Finally using this relation and the following Wronskian relation (page 489 12]), 
(Ifctffc-i + 4-iK fe ) (10 (*)) ' 



x (s) ' 
we obtain that 

2 



0D4 (x+ (s)) 
thus 

^ (a) , *) = nnrW 2 rrr 4 (a; (s) A x («)) (4 (s+ (*)) *4 - K (*+ (a)) h) (x (s) V x Q («)) . 
QDI k (x+ (s)) V J 
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We can now express the solution p for 8 < 9 by 



., I k (x (s) A x (s)) (l' k (x+ (»)) K k - K k (x+ (s)) I k ) (x (s) V x (s)) 

P 0, 0, s) = — > -7- — sin (k9) sin (k8 ) . 

GD ^ I k (x+ (s)) 



The exit point distribution e° (r) is given by 

t° (r) = - (j ^ (^°° P (r, (?, t) ) (0 = 0) . (3.20) 

To obtain an analytical expression for expression (|3.20p . we use the Laplace relation: 

F(z) 



CU f(u)du 

where F = L (/) is the Laplace transform of the function /. We have 
p(r,9,u)du = £- 1 fP^ r ' e ^y 



2 Ik 



k (x(s)Axo(s))(l' k (x + (s))K k -K' k (x + (s))I k )(x(s)\/xo(s)) 

- l- 1 [m5 E sin ^ sin ^ " 

The computation of the integral 

1 „ r+ioo I k (x(s)Axo(s))(l' k (x + (s))K k -K' k (x + (s))Ik)(x(s)y x (s)) 

I (r, 0, t) = — — ]T ™ (kO) sin (k9 ) / ^ e st d 

®kDi J_ too sl k (x+ (s)J 

(3-21) 

uses the residue theorem and the details are given in the Appendix. We have 
/ (r, 6,t) = Jp (r, 9, u)du=-^ (St (r, 9, t) + S 2 (r, 9,t)) , 

where 



Si(r, 9, t) = sin (k9) sin (k9o) 



2kR 2k 



'o 

OO 



o < a i\ ■ (ua\ ■ fua \V -Da 2 ,t J k (rotj t k) Jk (ro<Xj,k) 

b2{r, 9, t) = — 2 y sin (k9) sin (kvo) > e 



and Jfc are the fc-order Bessel's function and a^ft are the roots of the equation: 

J'k = o. 

Consequently, for r < ro, using (|3.20|) . we get the following exit distribution (for 
6 = 0): 



w = C lim e ' *) + ^ e - 

B raw Vt^oo 
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Because : 



lim Si(r, 9, t) = S^r, 9) and lim 5 2 (r, 0, i) = 0, 

t — >oo t — >oo 



we finally obtain that 

i „k— 1 ( r 2k I p2k\ 

e ° W = Q E ™ (**>) S (3-22) 

and, for r > ro, a similar computation leads to : 

i fe / 2fe , r>2k\ 

These expressions can be further simplified. Indeed, we rewrite them as follows (for 
r < ro) : 

e°(r) = ^sin(k9 Q ) fr ^^ ' 

thus, 



9r V Uo/ V Vi? 



where 3w denotes the imaginary part of the expression. We obtain two geometrical 
series that can be summed. We get: 



e° (r) = —3m 
Br 







git/So ( X. 


)" 










1 - e we ° ( 









1 - e 



00 feXjiT 



that is: 



6r 



1 - e l 



After some rearrangements, we obtain the following exit point distribution on 9 = 0, 
conditioned on the initial position (ro, 0o) : 

if (rr ) v sin (v6 ) {rr R 2 ) v sin(v9 ) 

e [r) = e u (r|r o ,0 o ) - ' 



6r yr 2 » + r 2v - 2 (rr f cos (i/0 o ) (rr ) 2,y + R Av - 2 (rr i? 2 ) iy cos (i/0„) 

(3.24) 

for < r < i?. Similarly, for 9 = 0, we obtain 

e/N o / I «s if (rr )" sin{v9 Q ) (rr R 2 ) v sin{v6 ) 

e u (r) = e w (r|r o ,0 o ) = 1 



6r ^r a " +rg" + 2 (rr ) v cob (i^o) (rr ) 2l/ + i? 4 " + 2 (rr i? 2 ) iy cos (i/0 o ) / 

(3.25) 

We notice that letting R tends to oo, we recover the expressions computed in the open 
wedge case f (!3.15p and (|3.16[1 ). 



COMPUTATION OF THE MEAN FIELD INDUCED BY MICROTUBULES IN CELLS 13 



£?-l ■■E'-ui-jnl'h 



in 
in 




Fig. 3.1. Mean exit points distribution. The theoretical distribution (dashed line) is 
tested against the empirical one (solid line) obtained by running a simulation of 20 000 Brownian 



particles, starting on the wedge bisectrix 



at ro 



R 



100 for 8 = ^ ), Because the 



starting point is located on the bisectrix, e° (x) = e e (x) , and thus the analytical curve is given by 

£ 8 („\ _ 2 ( (rr )("> , (rrpfl 2 )'"' 



2 
Or 



e (r) = e° (r) + e @ (r) : 
e (r) is achieved at r = roe 2v ' \ L '+ 1 ) . 



In t/ia4 case, t/ie maximum of the function 



3.4. The Mean Exit Radius (MER). To determine the mean exit distribu- 
tion radius e(r|ro) for a viral particle starting initially at position ro,#o where #0 
is uniformly distributed between and O, we consider e (r\r o,0q) = e° (?i?"o> $o) + 
e e (r\ro,0o) and estimate the integral 



e(r|r ) 



e 



e =o 



e(r|r , 60) Re- 



integrating expressions ( (13.241) and ()3.25|) ) we get 



(3.26) 



e(r|r ) 



r"-rff| 



+ Zn 



2u 



(rr ) 1 



i? 2,y - (rr )" 



We define the mean exit point as r m (ro) = E (r|ro) conditioned on the initial radius 
r . Thus, 



r m {ro) = E (r|r ) = / re(r\r )dr. 
Jo 



(3.27) 
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Using the expansion In (1 + a;) = J2 n >i (~ — for x < 1, we obtain by a direct 
integration that 



(2n+l)f „ 

e 



^,(2n+l) 2 (2 „ +1) \ (§) ^ ) j \ v ^ (2n+l)(((2n + l)§) 2 -l) 

(3.28) 



using the expansion in the first part, 

(2n+l) 2 (§) P=0 

and the approximation O << 1, we obtain using the value of the Riemann C— function, 
C(2) = £ and C(4) = f±, r < i?, that 



(r )«r 1 + — -— - (3.30) 



12 y Vi?y (Tr/ef-i 



For 8 small, the second term in the right-hand side of (|3.30|) is exponentially small. 

4. Approximation of a virus motion by an effective Markovian stochas- 
tic equation. We replace the successive steps of viral dynamics with an effective 
stochastic equation containing a constant steady state drift. 

4.1. Methodology. Virus motion described in paragraph (|2.2[) consists of a 
succession of drift and diffusing periods. We start with the stochastic equation 

X = —B^-r + v / 2~Dw, (4.1) 
|r| 

where r is the radial component of X , B is the amplitude of the drift. The MFPT 
of the process (|4.ip to the nucleus located r — S, when the initial position is located 
on the cell surface r = R is solution of 

K^ + ^) M, - s s (r - s, =- iforM) ^ 

t (r, (9) = for r = 8 
(r, 6>) = for r = R. 

dr 

A similar equation can be written in the domain Q with reflective boundary conditions 
of the wedge. Both processes in the full domain or in f2 lead to the same MFPT. The 
solution t(B,r) is given by 

j-R ( j-R -a(u-v) \ 

t(B,r) = C- / -^duU (4.2) 



where a = || and 



r R I pR a (u-v) \ 

t(B,R) = C = j I / —jyy du\dv. (4.3) 
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For a fixed radius R, the derivative of the function t (B, R) with respect to B is strictly 
negative, which shows that B — > t (B, R) is strictly decreasing. To determine the value 
of the amplitude B, we equal the mean time t (B, R) with the MFPT to reach the 
nucleus within the iterative procedure as described in paragraph (|2.2p : at time zero, 
the virus starts at a position r = R = Rq and reaches the edge boundary in a mean 
time u (Rq) and at a mean position r m (Rq). The viral particle is then transported 
toward the nucleus over a distance d m during a time t m . Either the particle reaches 
the nucleus before time t m and then the algorithm is terminated or in a second step, 
it starts at a position R\ — r m (Rq) — d m . The process iterates until the particle 
reaches the nucleus. We consider the mean number of fundamental steps (diffusion 
step and directed motion along a MT step) the virus needs to reach the nucleus is 
equal to n > 0. The mean time to reach the nucleus computed by equation (|4.2[) has 
thus to be equal to the mean time r = X)fc=o w(iik) + nt m + < t r > of the iterative 
trajectory. In a first approximation, wc neglect the mean residual time < t r > and 
we thus get the equality: 

n-l 

t(B,R) = T = Y J u(R k )+nt m (4.4) 

fe=0 

Rk+i = r m (Rk) - d rn (4.5) 
Rq = R . (4.6) 

For a fix radius R, equation (|4.4|) has a unique solution B, which can be found in 
practice by any standard numerical method. 

Remark. The MFPT of a particle where the trajectory consists of alternating 
drift (traveling along microtubules) and diffusion periods can either be higher or lower 
than the MFPT of a pure Brownian particle. Indeed when B < 0, the drift effect is less 
efficient than pure diffusion. For example, for Q = ^ , R — 100/xm, S = ^ = 25/xm, a 
large diffusion constant D = lOfim 2 s _1 with the dynamical parameters t m — Is and 
d m = 1/im, leads to a negative mean drift 

B w -0.14/OTis" 1 . (4.7) 

On the other hand, for a small diffusion constant D — l/im 2 s _1 , an efficient mi- 
crotubules transport obtained for t rn — Is and d m = 5/zm leads to a mean positive 
drift 

B w 0.13/xms" 1 . (4.8) 

4.2. Explicit expression of the drift in the limit of << 1. When the 
number of microtubules is large enough, the condition << 1 is satisfied. Moreover, 
because a virus entering a cell surface has a deterministic motion, we can assume that 
the initial position satisfies Rq < R so that we can neglect any boundary effects and 
use the open wedge approximation which consists of using formula (|3.30|) without the 
boundary layer term. Actually, this approximation is not that restrictive because after 
the first iteration process (movement along the microtubule followed by the particle 
release), the boundary layer term is negligible compared to the other term. 

To obtain an explicit expression for the amplitude B, we consider the successive 
approximations 

r m (R ) « ifo f 1 + fLV (4.9) 
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and 



Ro — Rq- 



e 2 



Ri ~ Rq [ 1 + — ) - d n 



02 

^- /! ' i ! 1 + T2 



e 2 

1 + l2 



fij — -Ro 1 



12 



£ 1 

\fc=0 v 



e 2 

12 



that is 



-Ri — hRo 



12rf„ 

e 2 



i 



e 2 

12 



e 2 



(4.10) 



Thus the particle reaches the nucleus after n iteration steps which approximatively 
satisfies R n — 6, 



h 



,ie- 

12d n 



i2o-<5 



Zn (1 + f£) d ™ 



o(l). 



(4.11) 



If T n denotes the mean time a viral particle takes to reach the nucleus, then using 
formula (13.911. we obtain 



^ tat 



(e) 



AD 



1 «-l 



2 



(4.12) 



that is 

n-l 



tan(e) 

e 



4L> 



E 

j=0 



12rf„ 

e 2 



12d„ ; 

"e 2 " 



i?0 — 



12d„ 

-Q2- 



»2\ ' 



e_ 
12 



i?n — 



12d„ 

"e 2 " 



12 



2 -A 



V 



/ tan{B) _ i \ 



AD 



(I2d n 



\ e 2 



24d r , 
-92 



-Ro 



12d„ 

-02 



1-1 



For 9 << 1, a Taylor expansion gives that 



12 



12 



-Ro_5 

dm 

(Ro - S) 

72D 



t m (Rq - 6 ) , #0 + (5 
24d m 



9 2 



d m + 3 (i? + 5) + 



2 (i? 2 + i? <5 + 5 2 ) 



Ro 



124 
-92 



2 1 



i-d+n; 



e 4 + o (e 4 ) . 
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In small diffusion limit D << 1,0 << 1, the velocity is B ~ s and consequently 
we obtain for Rq m R, a second order approximation 

Bk 7 ^fr^ > ( 4 - 13 ) 



l 



where d m ,t rn are the mean distance and the mean time a virus stays on the micro- 
tubule, R (resp. 6) is the radius of the cell (resp. nucleus) and = tS, where N is 
the total number of microtubules. 

4.3. Justification of the MFPT-criteria.. To justify the use of the MFPT- 
criteria to estimate the steady state drift, we run numerical simulations of 1,000 
viruses inside a two dimensional domain Q (6 < r < R) with intermittent dynamics, 
alternating between epochs of free diffusion and directed motion along microtubules 
and compare the steady state distribution with the one obtained by solving the Fokker- 
Planck equation for viruses whose trajectories are described by the effective stochastic 
equation (|2.2|) with our computed constant drift 

d m 

b(X) = r -L = -.B-1 (4.14) 

We imposed reflecting boundary conditions at the nuclear and the external membrane. 
The theoretical normalized steady state distribution p satisfies 

DAp — V.[bp] = in Q, 

and the solution p is given by 

Br _ Br_ 

D ED 

p(r) = _ — = , (4.15) 

J s e d 2irrdr 2vr§ ( 5e-~ - Re-— + § ' - ~~ 



The result of both distributions is presented in figure 14.11 where we can observe that 
both curves match very nicely. This result shows that the criteria we have used is 
at least enough to recover the distribution. For the simulations, we consider the 
directed run of the virus along a MT (loaded by dynein) lasts t m = Is and covers 
a mean distance d m — 0.7p,m 13J. The diffusion constant is D = 1.3pm 2 s^ 1 as 
observed for the Adeno Associated Virus [14]. The two curves in figure [4~T1 fit very 
nicely except at the neighborhood of the nuclear membrane, where the simulation of 
the empirical distribution is plagued with a possible boundary layer. Another source 
of discrepancy comes from the difference of behavior of viruses far and close to the 
nucleus: viruses far from the nucleus do not bind as often as those located in its 
neighborhood. Consequently, a constant effective drift cannot account for the radial 
geometry near the nucleus. A theory for radius dependent effective drift has been 
derived in [15] . 

5. Conclusion. In the limit of a cell containing an excess of microtubules, we 
have presented here a model to describe the motion of biological particles such as 
viruses, vesicles and many others moving inside the cell cytoplasm by a complex com- 
bination of Brownian motion and deterministic drift. Our procedure consists mainly 
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Fig. 4.1. Steady State distributions. We show the empirical steady state distribution for 
1,000 viral trajectories with an intermittent dynamic (solid line). The theoretical distribution of 
viruses whose trajectories are described by the stochastic equation (2. 2V is given in dashed line. 
Geometrical parameters are : R = 20fim, 8 = 5jim and © = t?z . 



in approximating an alternative switching mode between diffusion and deterministic 
drift epochs by a steady state stochastic equation. This procedure consists of estimat- 
ing the amplitude of the effective drift and is based on the criteria that the MFPTs to 
the nucleus, computed in both cases are equaled. In that case, this amplitude account 
for the directed transport along microtubules, the cell geometry and the binding con- 
stants. The model has however several limitations. First, we do not take into account 
directly the backward movement of the virus along the microtubules [TH1 [T7] , which 
can affect the mean time and the amplitude of the drift. Second, the present com- 
putations are given for two dimensional cell geometry only. It can still be applied to 
many in vitro culture cells, however it is not clear how to generalize our approach to 
a three dimensional cell geometry. For example, to study the trafficking inside cylin- 
drical axons or dendrites of neuronal cells, a different approach should include this 
geometrical features. However despite these real difficulties, the present model may 
be used to analyze plasmid transport in an host cell, at the molecular level, which is 
one of the fundamental limitation of gene delivery [T51 UM HZT] . 

Appendix. In this appendix, we provide an explicit computation of integral 
(|3.2ip using the method of the residues. This method was previously used in a similar 
context in ([T^] p 386). We denote by (Pj)j >0 the poles of the function 

I k (x (s) A x (a)) (l' k (x+ (s)) K k - K k (ar + («)) I k ) (x (s) V x (s)) 

$ : s _> i . L e s *. 

sI 'k ( x + ( s )) 

where (x (s) = r^p^ , xq (s) = r y^ and x + (s) = Ry^). The associated residues 
are (r^) , >Q . We now compute the residues explicitly. 
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To identify the poles, we recall the relation between the fc-order Bessel's function 
J k (that is true for z such that — 7r < arg (z) < ^) and the modified Bessel functions 
4 (P 375 [Til): 

I k {z) = e- L - km J k (ze^y (5.1) 
All roots a^k of the equations 

J'k (Ra) = 0, 

are real, simple and strictly positive (p 370 [TT]) because k is real and 

k < ax,k < oi 2 ,k ■ ■ ■ 

Thus, 

l'k (- iRa 3,k) = 0. 

Finally the poles of $ are simple given by p§ = and Vj > 1, p| = —D<Xj k - Conse- 
quently the associated residues are given for each k for all j > by 

rj = lim (s-p)) $(a). (5.2) 
Then using the residues, integral (|3.21[) is given by 

1 (r ' $ ' t] = OrrDi ^ Sm m ( ^ o) (2 ™ } ^ = 4d ^ {W) ^ ( ^ o) ^ ^ ' 

We now compute the residues . The residue r§ is associated with the pole = 
and given by 

7-q = lim s<&(s) 

s— »0 

Using the following identities on the modified Bessel functions (p 489 [T2]) 

1 k 1 k 

4 (z) = 4+1 (z) + -I k (z) and K k (z) = -K k -i (z) - -K k (z) , 

z z 

substituting the derivatives I k and K k in the expression of <&, we get 

k v h(x{s) AzoO)) 

= hm 



(4+1 + s^jlfc) (z+ 00) 

Taking into account the dominant terms only, we get 

fe _ ,. 4 (x (5) A x (a)) (4 (g+ (g)) gfc + J4 (g+ (s)) 4) (a; (5) V z (3)) 
—0 4 ( x+ («)) 
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To further compute this limit, we use the Taylor expansions of I k and K k (p 375 [11) ) 
expressed in terms of the T function: 

J* W « and Jir fcW «ir (*)(!,) \ 

For r < ro, we get 

%S# (%# *r (*) (| (xo (.))) fe + |r (*) (i {x+ is))) k 

>"n = lim ; u 

s^O (%( x+ ( s ))) k 

r(fc+i) 

Finally, using the relation r (k + 1) = fcT (fc), and the expressions of x(s), Xq(s) and 
x + (s) we get 

fc _ r k (rf + R 2k ) 
r ° = 2kR 2k r k ■ 

The computation of the other residues { r j)j >1 > is slightly different 

r k = lim (s-Pj) $(s), 

where p k — —Da 2 k . Using the Wronskian relation (p 489 [12]) : 

I k {z)K' k {z)-K k {z)I' k {z) = -- 

z 

we now substitute 

, -l + K k (z)l' k ( Z ) 

Kk[z) - hJJ) ' 

in the expression of $, we get 



lim 

; s I k (x+{s)) 



(s _ pk) &st h (x (*)) \I h (x + (s)) K k - -j (x + («)) I k j (x (»)) 



Because 



lim l' k (x + ( S )) = l' k (x + (p k )) =0, 



we obtain the expression for the residues: 



^ h(x(p k ))h{x,{p k )) ^ (s-p?) 
Pi 4- (x+ (pf)) z+ (p£) s^v) I' k {x+ (s)) ' 



Finally, since 



( s _ pfc) _ 2^ Dp) x+ (s) _ ^ ( p fc) 2^; 



lim — = lim 

s ^l' k (x + (s)) R s^ P > l' k (x + ( S ))-l' k (x + (p k )) Rl' k (x + (p)))' 
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we obtain 

e# !*(*(,,*))/* (so (#)) 2 ^ 



r k 

' 3 



To simplify this expression, we use that Ik satisfies the differential equation (p 374 

mi): 



z \ z 



thus for z — x + (p 1 ?) 



ii (* + w)) = pf ^f fc2 4 (* + m . 

we get 

fe _ 2De& I k (x(tf))l k (xo(rf)) 



i? 2 ^ + Dk 2 


3 (*+(#)) 


get 








-i? 2 4 fe + fc 2 


Jl {Ra 3 ,k) 



r k = 

' 3 

Integral (|3.2ip is given by 

I(r, 9, t) = Yl sm ^ sm E r " = -§D {Sl{r ' °' t] + 5a(r ' <)} ' (5 - 3) 

k j>0 

where 

r fe ( r 2fe + R 2kj 



Si(r, 9, t) = sin {k9) sin (k6o) 



2kR 2 



k,,,k 



•0 

OO 



S 2 (r, 9, t) = -2^2 sin ( ke ) sin ( fc6 W E 



3 = 1 



(R 2 al k -k 2 )jl{Ra hk ) 
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